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Abstract 

A computational scheme is described which is second-order accurate m 
time and fourth-order accurate in space (2-4). This method is applied to 
study the stability of compressible boundary layers. The laminar compressible 
Navier-Stokes equations are solved with a time harmonic inflow superimposed on 
the steady state solution. This results in spatially unstable modes. It is 
shown that the second-order methods are inefficient for calculating the growth 
rates and phases of the unstable modes. In contrast the fourth-order method 
yields accurate results on relatively course meshes. 
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I . Int roduc t ion 


This paper is concerned with a fourth-order accurate finite difference 
scheme for the compressible, unsteady Navier-Stokes equations. The primary 
interest is the computation of spatially unstable disturbances into the 
nonlinear regime and the active control of such disturbances. Due to the 
wavelike nature of these disturbances, this problem has features of both wave 
propagation and fluid dynamics and the numerical scheme must be chosen to 
accurately compute waves propagating in an unstable, high Reynolds number mean 
flow. The application of this scheme to study the active control of spatially 
growing disturbances has been described previously [1]. This paper describes 
the numerical scheme and the advantages that can be obtained by the use of 
fourth-order accuracy. 

Higher order accurate methods, in particular spectral methods, have been 
successfully used in the computation of incompressible flows. Examples of 
such calculations can be found in [2] - [4]. Generally, spectral methods 

assume that problem is periodic in the streamwise direction so that Fourier 
(as opposed to Chebyshev) collocation can be used. (Spatially unstable 
disturbances were considered, however in [4].) The use of higher order 
methods for the numerical computation of waves has also been extensively 
considered in the literature. Examples can be found in [5] - [7]. 

In section 2, we describe the numerical scheme and discuss certain 
implementation details. In section 3 the performance of this scheme on a 
variety of problems is illustrated. Finally, in section 4 we discuss these 


results and draw conclusions. 
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II* Numerical Scheme 

The numerical scheme is an extension of the second-order MacCormack 
scheme due to Gottlieb and Turkel [8]. For the one-dimensional equation 
U t = F^ where F = F(U), we have (letting the superscripts denote the time 
level and the subscripts denote the spatial grid points) 
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(2.1a) 
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There is an obvious symmetric variant of (2.1) starting with a backwards 

predictor and then a forwards corrector. 

If F depends only on U it is shown in [8] that the scheme (2.1) is 

second-order accurate in time and fourth-order accurate in space ((2-4) 

scheme). Specifically, if At is the time step and Ax the space step, then the 

4 2 

truncation error is 0(At((Ax) + (At) )J. Thus the scheme is fourth-order 
accurate provided At = 0((Ax)^) as Ax -► 0. For nonlinear problems this is true 
only if the two variants of (2.1) are alternated. 

Although true fourth-order accuracy is obtained only for At = 0((Ax)^) it 
has been found that (2.1) is considerably more efficient than second-order 
schemes (see for example [6], [9-10]). For two-dimensional problems (2.1) can 
be used together with operator splitting [11] to maintain the (2-4) accuracy. 
Specifically, for the equation 


U - F + G 
t x y 


we have 
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U n+2 = L L L L U n 
x y y x 

where L x and Ly are one-dimensional solution operators corresponding to the 
scheme (2.1) applied to the equations U t = F x and U t = Gy respectively. 
Alternatively, one can consider an unsplit version of (2.1) [12]. 

The explicit nature of the scheme (2.1) makes it well suited for current 
vector computers such as the CDC Cyber 205 and the Cray 1-S and XMP series. 
The scheme has been implemented on the CDC VPS 32 at NASA Langley Research 
Center using vector operations over a two-dimensional grid (vector lengths 
20,000). The explicit nature of the scheme makes it relatively inefficient 
to compute steady flow, unless acceleration techniques are incorporated 
[13]. For the unsteady flows considered here, a time step restriction is 
necessary in order to accurately resolve the fluctuations in time and so 
explicit schemes become efficient. 

The implementation of the scheme (2.1) is straightforward if the flux 
function F depends only on the unknown U. For the Navier-Stokes equations we 
have F = F(U,U x> Uy) and in the evaluation of F^ it is necessary to approximate 

U x [8]. For a forward sweep, i.e., (2.1a) U x is approximated by a two-point 

backwards difference, i.e., U x ~ - U^_^/Ax and conversely for a backwards 

sweep. It is shown in [8] that for a general, non-constant coefficient 

problem this results in a scheme that is only third-order accurate. However, 
it can be shown that the third-order truncation error is eliminated by 
alternating the sweeps and the resulting scheme is fourth-order accurate for 
both the inviscid and viscous terms. Mixed derivatives, i.e., terms of the 
form Uy in an x sweep are approximated by second-order central differences. 
In many applications these terras are relatively small. We have not 
experimented with fourth-order differences for these terms. 
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In order to complete the description of the numerical scheme, it is 
necessary to describe the implementation of boundary conditions. We 
distinguish between boundary conditions which must be imposed as part of the 
problem and those boundary conditions which are necessary because the 
straightforward application of the difference formula (2,1) is not valid at 
boundary points* In this section we consider only the latter. 

Consider the forward predictor (2.1a). If l = N denotes a boundary point 
then the values of are not available for i > N and thus the scheme (2.1) 
can not be applied at the points i = N-l and i = N. The most satisfactory 
boundary treatment that we have found is to define F x for i > N by third-order 
extrapolation from the interior (see [6] and [8]). Hence, and F^ + 2 are 
defined by 


N+l 


= 4F - 6F 

N-l 


+ 4F 


N-2 


- F. 


N-3 


( 2 . 2 ) 
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= 4F 
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In (2.2) the extrapolated fluxes include both the viscous and inviscid terms. 

The scheme (2.1) has a five-point stencil. Implicit fourth-order schemes 
with a three-point stencil can be constructed by using Pade approximations. 
It is well known that compact schemes are more accurate than five-point 
schemes [14], however, they are more expensive because of the implicit nature 
of the scheme. One way to increase the accuracy of the scheme (2.1) is to use 
a scheme which is sixth-order accurate in space. The scheme (2.1) can be 
easily extended to sixth-order. The resulting scheme is 



37 f F n - F n l - 8fF n - F n ) + fF n - F n ) 
^ l+l 1 i+2 i+l^ 1 i+3 i+2 J 
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The suitability of (2.3) for the unsteady Navier-Stokes equations is currently 
being investigated. 


Ill • Numerical Examples 

In this section we describe some numerical examples illustrating the 
effects that can be expected from the fourth-order accurate differencing. The 
present program is primarily designed to compute unsteady flows and does not 
make use of any acceleration techniques which destory the consistency in time 
[13]. The first numerical example will, however, illustrate the effect of the 
higher-order differencing on steady flows. 

Specifically, we consider a supersonic boundary layer over a flat 
plate. The free stream Mach number is 4.5 and the unit Reynolds number is 
2.3 x 10^. The inflow data is generated by a boundary layer program [15] at 
0.5 ft. from the leading edge. The computational domain is 2.0 ft. in x 
and 11 6 in y where 6 is the boundary thickness at inflow. An exponentially 
stretched grid is used normal to the plate with the Jacobians evaluated 
analytically. Figure 1 illustrates the computational domain. 

In addition to the numerical boundary treatment described previously, it 
is necessary to discuss the boundary conditions which must be imposed. For 
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the supersonic case all quantities are imposed at the inflow while the 
solution at the outflow boundary is obtained by zeroth-order extrapolation 
from the interior. We have verified that the solution is very insensitive to 
the treatment of the outflow boundary as would be expected since the 
linearized characteristic variables are convected out of the computational 
domain. At the plate, the two velocity components are set to zero and the 
temperature is specified. Several different techniques have been used to 
obtain the pressure. These include different orders of extrapolation in the 
normal direction and the use of the normal momentum equation. The accuracy of 
the solution is insensitive to the boundary treatment (probably because of the 
stretched grid). At present we simply use first-order extrapolation for the 
pressure. 

The treatment of the upper boundary (see Fig. 1) is based on the 
linearized characteristics m the normal direction. 

The Navier-Stokes equations can be written in the generic form 


U = F + G 
t x y 


(3.1) 


T 

where U is the vector (p, pu, pv,E) , p is the density, u and v are the x and 
y velocities, respectively and E is the total energy. The functional forms of 
F and G are standard and are omitted for brevity. In the normal direction we 
consider only the system 


U 


t 



(3.2) 


In the applications the vertical velocity v is small and positive at the upper 
boundary. Upon linearizing (3.2) we find that the three characteristic 
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(3.3a) 
(3.3b) 
(3.3c) 

are convected to the boundary from the interior. These three variables are 
obtained by zeroth-order extrapolation from the interior. Since we consider 
the linearized characteristics, the quantities with a " in (3.3), are taken 
from the previous time step. The final boundary condition is obtained from 
setting 

P t ~ ^ V t = ° (3.3d) 

corresponding to the characteristic variable entering the computational 

region* The use of the radiation condition (3* 3d) was previously found to 

considerably accelerate the convergence to the steady state and to permit the 

upper boundary to be brought closer to the plate [16]. 

A typical comparison between the second and fourth-order schemes is shown 

in Fig. 2. The streamwise velocity is plotted against y at a location of 1.0 

ft. from the leading edge (Re * * 17174), where Re . is the Reynolds number 

6 5* 

based on displacement thickness. The results are shown for the fourth-order 
scheme (2.1), the second-order MacCormack scheme and the solution obtained 
from the boundary layer equations. The finite difference results were 
obtained from the same grid (21 x 31) and the improvement in accuracy of the 
fourth-order scheme is evident. It should be noted that we do not obtain a 
similar improvement for more viscous boundary layers. The reason for this is 
not understood and is currently being investigated. 


variables 


p - pc 
p + pcv 
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For our next example, we consider an unsteady disturbance in a relatively 

low speed subsonic boundary layer* The mean flow is a boundary layer with 

free stream Mach number 0*4 and a unit Reynolds number of 3.0 x 10^. At 

inflow we have Re * = 998 and the computational domain is chosen so that at 
6 

outflow Re * = 1730. A fluctuating disturbance is introduced at the inflow 
6 

with nondimensional frequency F = ( 27 ifv)/u£ of .8 x 10^. f is the frequency 

m Hertz, v the kinematic viscosity, and the free stream velocity. This 

flow is nearly incompressible. Based on linear, incompressible stability 

theory it is known that this frequency is unstable at inflow but becomes 

stable at Re * = 1450. Since this is a subsonic flow, three boundary 
6 

conditions must be imposed at the inflow- These conditions are obtained by 
computing an eigenfunction of the Orr-Sommerf eld equations (using a program 
developed by R. Dagenhart of NASA Langley Research Center). The real part of 
the solution (times e^ Ft ), is then used to compute the three linearized 
characteristic variables which enter the computational domain from the 
outside • Specifically, 


p + pcu 
P - pc 2 
v 

where e determines the strength of the disturbance. The outgoing 

characteristic variable p - pcv is obtained from the interior by zeroth-order 
extrapolation. 

The outflow boundary condition is treated similarly. The incoming 
characteristic variable is set equal to zero by imposing the condition 


= mean + e (values obtained from Orr-Sommerf eld equation) (3.4a) 
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P t - pcu t = 0, (3.5) 

The use of the condition (3.5) and more accurate radiation condition to 
compute unsteady disturbances in subsonic flows is discussed in [16]. 
Characteristic radiation boundary conditions of the form (3.5) are commonly 
used m two-dimensional linear wave propagation problems (for example, see 
[17]). The boundary condition (3.5) has been found to be sufficiently 
accurate for the present calculations. This is based on comparisons with 
linear stability theory. 

In Fig. 3 we plot the computed growth rate as a function of Re *. The 

6 

/ 2 2 

growth rates are computed by calculating the RMS =y (pu) - (pu)^ ^ anc * 
integrating in y. The inflow data is chosen so that the maximum perturbation 
is 2% of the free stream. In this case the (2-4) scheme shows a significant 
amount of nonlinear growth while the (2-2) scheme, for the same grid, shows 
significantly less growth. Mesh refinements verify the accuracy of the fourth- 
order code. 

Reducing the inflow perturbation to 0.2% of the free stream the fourth- 

order results reproduce the linear theory results very closely and in 

particular, show that the perturbation becomes stable at Re * 1450 as 

6 

predicted by linear theory. 

In Figs. 4a and 4b we plot pu as a function of time for Re * = 1263 and 

6 

y = 0.0034 ft. In Fig. 4a the (2-4) scheme is shown for two different grids 

and shows that a further grid refinement does not yield any additional 
information. In Fig. 4b the (2-2) scheme is compared with the (2-4) scheme. 
At this location the disturbance is basically linear and the figures indicate 
significant amplitude and phase errors for the (2-2) scheme. 
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Figures 5a and 5b contain similar results further downstream at the 

location Re * = 1481 and y = 0.0011 ft. At this location the disturbance 

6 

exhibits some nonlinear effects which are not found by the second-order scheme 
with this mesh. 

In Fig. 6 the computed growth rates are shown for a disturbance m higher 

speed flow where compressibility effects would be expected to be important. 

The free stream Mach number is 0.7 and the unit Reynolds number is 300,000. 

The computational domain is chosen so that at inflow we have Re * = 900. The 

6 

nondiraensional frequency F = 0.8 x 10^. The inflow data was generated from 

the incompressible stability program, however, we have not compared growth 

rates with those predicted by the incompressible stability theory. The 

results in Fig. 6 show a significant reduction in the growth rate using the 

second-order scheme. Finally, in Fig. 7 we plot the mean of pu across the 

boundary layer for Re * = 1400. It is apparent that the second-order scheme 

6 

misses significant features of the flow. 


V • Conclus ion 

There are several types of errors that appear in the computation of 
unsteady waves and stability analysis. These include both amplitude and phase 
errors. These errors frequently simulate an apparent lower Reynolds number 
and indicate that the flow is more stable than it physically is. These 
effects are seen in the results presented here. It is shown that second-order, 
m space, schemes are not efficient m simulating nonlinear stability of high 
Reynolds number flows. The existing program solves the two-dimensional 
laminar compressible Navier-Stokes equations and is currently being extended 
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to three dimensions. Both subsonic and transonic stability regions can be 
calculated using the improved fourth-order method. 
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Fig. 1. Schematic of computational domain. 
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Fig. 3. Comparison of amplitude growth for the second 

and the fourth-order schemes. M =0.4. 

00 



Fig. 4a. pu vs. time for two grids at Re 


= 1263, y = 0.0034 ft 
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Fig. 4b. pu vs. time using second and fourth-order schemes. 

Re *= 1263, y = 0.0034 ft. 
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Fig. 6. Comparison of amplitude growth for the second 
and the fourth-order schemes. M =0.7. 
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